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Abstract 

We present a method for performing time domain simulations of a microphotonic system con- 
taining a four level gain medium based on the finite element method. This method includes an 
approximation that involves expanding the pump and probe electromagnetic fields around their 
respective carrier frequencies, providing a dramatic speedup of the time evolution. Finally, we 
present a two dimensional example of this model, simulating a cylindrical spaser array consisting 
of a four level gain medium inside of a metal shell. 
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I. INTRODUCTION 



Interest in microphotonic lasing systems has been increasing over the past few years. As 
a result, it has become more important to be able to numerically simulate these lasing sys- 
tems. Several finite difference time domain (FDTD) simulations of a four level gain medium 
embedded in a microphotonic system have been presented previously^"-, but these simula- 
tions all use structured (cubic) grids and consequently accurately model curved geometries. 
There have also been methods developed to model spherical gain geometries by expanding 
electromagnetic fields as sums of spherical Bessel functions^^. These methods overcome 
the limitations of structured grids for spherical geometries, but in turn are limited to only 
modelling spherical geometries. In principle it should be possible to model a microphotonic 
lasing system with an FDTD simulation utilizing unstructured grids, but to the best knowl- 
edge of the authors this has not been demonstrated. The finite element method (FEM) can 
utilize unstructured grids and as a result can model a wide variety of geometries. In this 
paper we present a FEM model of a microphotonic system with gain arising from a four 
level quantum system. In addition to developing a FEM microphotonic lasing model, an 
approximation is introduced whereby the pump and probe fields are solved for separately, 
with each field described by the slowly varying complex valued field amplitude of a con- 
stant frequency carrier wave. This approximation allows for much larger time steps and a 
considerable speedup in simulation time. 

In the first section of this paper we will describe the dynamics of the microphotonic 
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lasing system, the carrier wave approximation, and finally the finite element formulation of 
the problem. In the second section we present a two dimensional model of a one dimensional 
cylindrical spaser array as an example of this new simulation method. 



II. FEM MICROPHOTONIC LASING SIMULATION 
A. Field equations of a microphotonic lasing system 

The simulation we present of a microphotonic lasing system requires the time domain 
modelling of several different fields and their mutual interactions. These fields include the 
electromagnetic field, the electric polarization field inside a metal with a Drude response, 
the electric polarization field of the gain medium, and the population density fields of the 
different energy levels of the gain medium. Each of these fields evolve according a particular 
differential equation that must be solved when simulating a microphotonic lasing system. 
The field equation for the electromagnetic field is 



where A is the electromagnetic vector potential, P is a polarization vector describing either 
a Drude response from a metal inclusion or a Lorentzian response from a four level gain 
system, and e r is a relative permittivity that is constant with respect to frequency and not 
included in P. Here and for the remainder of the paper we have use SI units. Also, we 
have used the temporal gauge condition dA /dt = along with the initial condition for the 
electrostatic potential A (t = 0,x) = to ensure that the electrostatic potential is zero for 
all time, eliminating it from our equations. Given this choice of gauge the electric field and 
magnetic flux density are defined as 



B = V x A 

Using this definition for E, the Drude response of a metal inclusion is determined by the 
equation 




(1) 
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+ 7 P = -e c^A 



(3) 



where u p is the plasma frequency and 7 is the damping frequency of the Drude metal. 
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FIG. 1: Simple model of a four level gain medium. The lasing and pump transitions are assumed 
to be electric dipole transistions with frequencies of oj a and uib respectively. The decay processes 
between the i-th and j-th energy levels are described by the decay rates l/Ty. 

The gain medium is modelled as simple four level quantum system, described schemat- 
ically in Fig. [U The 1—^2 transition is an electric dipole transition with a frequency of 
bj a . Similarly, the — > 3 transition is also an electric dipole transition with frequency Ub- 
Spontaneous decay between the i-th level to the j-th level occurs at the decay rate of l/r^-. 
These decay rates include both radiative (spontaneous photon emission) and non-radiative 
(spontaneous phonon emission) decay processes. In the case of spontaneous photon emission, 
our model does not produce a photon. Coupling of the gain medium to the electromagnetic 
field is only allowed for stimulated photon emission. 

The electromagnetic response of the four level gain system is given by 

^f + r ^+^ai = -* B (N* - N l4 )E 4 , 

^|# + r^ + ^ 2 P bj = -<r 6 (N 3i - NoOEj. 
Here P° and are the i-th components of the gain polarization due to transitions between 
the 1st and 2nd levels and between the Oth and 3rd levels respectively. Additionally, T a and 
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Tb are the linewidths of these transitions, a a and cr& are coupling constants, and N j, Nh, N 2 j 
and N 3 j are the population number densities for oscillators polarized in the i-th direction for 
the Oth, 1st, 2nd and 3rd energy levels. Note that T a > l/r 2 i and T b > l/r 3(f ^. 
Finally, the population number densities evolve according the equations^ 1 ^ 



(5) 



dt hu b dt \r 30 r 32 
dN^ = Nsi 1 E dP ai N 2l 

dt T 32 ku a dt Til ' 

9Nu = N 2i 1 E dP ai N u 
dt r 2 i huj a % dt no ' 

gN 0t = N 3i Ni, 1 E dP bl 
dt r 30 no ^o; 6 * <9t 
Together, this system of equations (Eqs. ( TTl[3l l5|)) completely describes the dynamics of the 

microphotonic lasing system. The main disadvantage of solving this system of differential 

equations is the small time step required. In practice, ~ 100 time steps per period of the 

pumping laser beam are required for an adequate simulation. A typical lasing simulation 

could require over 100,000 lasing periods, making the computational requirements of the 

simulation prohibitively large. 



B. Period averaged approximation 

There is a simple method for dramatically speeding up the simulation time. The elec- 
tromagnetic field as well as the polarization fields oscillates at two frequencies. These two 
frequencies are approximately equal to the frequency of the 1 — > 2 transition frequency u a , 
and the — > 3 transition frequency Ub- Much of the computational effort required in this 
time domain simulation is spent on these simple, approximately harmonic oscillations. A 
good approximation is to assume these fields oscillate harmonically, with complex valued 
amplitudes that are slowly changing in time. We can ignore the fast oscillations and instead 
simulate the relatively slower time dependence of these amplitudes. 

Since there are two frequencies, we divide our electromagnetic field into two separate 
fields 
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A(t, x) = A,(t.x)^' + A,( t ,x )e ^ + c.c._ (6) 

with each field oscillating at a different frequency. Here Ai is the complex valued amplitude 
for an electromagnetic field that oscillates at a frequency close to the 1 — > 2 transition 
(ux ~ (J a), and A2 is the complex valued amplitude for an electromagnetic field that oscillates 
close to the — > 3 transition (u; 2 ~ cu&). Also, c.c. indicates the complex conjugate of the 
preceding terms. 

By inserting the above equation into Eq. ([1]), the field equation for A, we derive two new 
field equations 

„ / 1 \ f q 0Ai 9 2 AA <9Pi 

V x —V x Ai + e r e -u(Ax + 2iwi— — + 

'' ct~ / at 

(7) 



/i V V dt dt 2 J dt 



V x ( —V x A 2 ) + e r e ( -wf A 2 + 2iw 2 ^ + 



v ' u v 1 dt dt 2 J dt 

Here we have also separated the polarization field into two fields 

p(t x) = Pi(t^)e i ^ t + P 2 (t^)e i ^ t + c.c (g) 
For Drude metal inclusions the polarization fields obey the equations 

iu 1 P d 1 + — i+ 7 P? = -e ^Ai, 



<9P d 

\uj 2 P d 2 + -^ 2 + 1 P d = -e u 2 p A 2 , 
while for Lorentzian gain inclusions the polarization fields obey the equations 



(9) 



f>p9 f)2p9 f Op9 

-u 2 P 9 2l + 2ic 2 ^|i + + T b (iu 2 P 9 2l + ) + "iKi = ( N* - No,- ) E 2/ , 



(10) 



Here and E 2i are the i-th components of the electric fields associated with the potentials 
Ai and A 2 , and are defined as Ei = —dAi/dt and E 2 = —dAi/dt respectively. 
Finally, the new differential equations for the occupation number densities are 
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(11) 



hub \ 1 dt 



r 30 r io 

Here the coupling term between the occupation number density fields and the electromag- 
netic and polarization fields has been replaced by a term representing the period averaged 
value of these terms 
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where * indicates complex conjugation. 

Finally, we mention that while we have developed the preceding approximation using a 
FEM model, this approximation is not limited to the FEM. It could potentially be used 
to speedup both the FDTD models of microphotonic lasing systems^- as well as the time 
domain models of spherical lasing geometries utilizing spherical Bessel functions^ 1 ^. 



C. Finite element formulation 



Now that we have derived the period averaged field equations (Eqs. (|7f9til2p ) for the 
microphotonic lasing system, we can convert these differential equations into weak forms 
that can be solved in a finite element simulation. 

The weak forms for the field equations of the electromagnetic fields (Eq. (J7])) are 



S 



Al 



;Ai, Ai) = (vxAi)-^(Vx Ai) + e r e Ai ■ (-u^Aj. + 2^1^ + 9 A 



dt 2 



" Al ' ~dT' 



(13) 
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A 2 , A 2 ) = (V x A 2 J • — (V x A 2 ) + e r e A 2 • I -u£ A 2 + 2iu 2 -^- + 



dt 2 



A dP2 

Here ~ indicates a test function^ 1 ^. These weak forms enforce both the electromagnetic 
field equations as well as a natural boundary condition^^. The finite element method 
requires that the integral of the weak form over the simulation domain be set to zero. As 
an example, if we apply this requirement to the weak form F^i, we find that by integrating 
by parts we obtain a volume integral enforcing the electromagnetic field equation as well as 
a second boundary integral enforcing a boundary condition on the field, 



(fx F 



Al 



d 3 x Ai • 



/'o 



V x —V x Ai + e r e -u^Ai + 2^—— + 



0Ai <9 2 Ai 
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dt 2 
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dA Ai 
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n x ( —V x Ai 



(14) 



Here Q is the simulation domain, dfl is the boundary of that domain, dA is a infinitesimal 
differential area on that boundary, and n is the direction normal to the boundary. In the 
absence of any extra boundary terms, the boundary integral in Eq. (114p forces the tangential 
component of the magnetic field Hi to zero. This perfect magnetic conductor boundary 
condition is not desirable for our simulation, so we will modify it to allow for a boundary 
that absorbs and emits plane waves at normal incidence to the boundary. 

For a flat boundary at a large enough distance from the inclusions in the simulation 
domain that evanescent waves are negligibly small, if the remaining propagating fields are 
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normal to this flat boundary then the vector potential can be represented as the sum of two 
vector potentials, 

Mt, x) = a - + b (t + ^) . (15) 

Here a is the vector potential of a plane wave propagating toward the boundary and b is 
the vector potential of a plane wave propagating away from the boundary. The boundary 
condition we desire is one that absorbs a and emits an arbitrarily defined b. If we take the 
part of the surface integrand from Eq. [TH that is within the brackets and substitute Eq. [15] 
for Ax we get 




where E^"* = —da/dt is the part of the electric field associated the plane wave a propagating 
toward the boundary, and E* 1 nc = —db/dt is the part of the electric field associated with the 
plane wave b propagating away from the boundary, the sum of which is E° u * + E* 1 nc = E x = 
—dAi/dt. Also, zo = a/a^oAo is the impedance of free space. Multiplying Eq. [16] by a test 
function Ai and integrating over the domain boundary gives us a new boundary weak term 

S A i(Ai,Ai) = - <f dA -Ax • 
J an z o 

Adding this additional boundary weak term to specific boundaries enforces a matched bound- 
ary condition (referred to as an absorbing boundary condition in Ref.— ) which allows for 
plane waves normal to the boundary to be absorbed and for the incident plane wave E™ c 
to be emitted into the domain normal to the boundary. A matched boundary condition for 
A2 can be enforced in the same manner. 

The weak forms for the remaining field equations are simpler since these differential 
equations only involve derivatives with respect to time. The weak form for the polarization 
of Drude metal inclusions is 
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where again a ~ indicates a test function. Similarly, the weak form for the polarization 



fields of the gain medium are 
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and the weak forms for the population density rate equations are 
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where the period averaged values for the coupling term are given in Eq. ( TT2|) . Also, we can 
avoid solving for N i by taking advantage of the fact that N 0i = N int — Ni» — N 2 j — N 3i where 
Nj ni is the initial value of No« when = N2i = N 3 j = 0. 

III. CYLINDRICAL SPASER ARRAY 



As an example of a microphotonic lasing system simulation we present a two dimensional 
model of a spaser (surface plasmon amplification by stimulated emission of radiation^^). 
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The time domain FEM simulation was performed using the commercial software COMSOL 
Multiphysics 3.5. For time stepping, the Generalized-a method was used with the damping 
parameter p in f = 1. A copy of the model can be obtained by contacting the corresponding 
author by email. 

The spaser is a one dimensional array of cylinders, each cylinder being infinite in extent 
in their axial direction. Each cylinder has a core consisting of a four level gain medium with 
a radius of r\ = 30nm and an outer shell composed of Ag with an outer radius of r 2 = 35nm. 
A diagram of the simulation domain is provided in Fig. |2j The artificial gain medium is 
characterized by the lifetimes r 10 = lCT 14 s, r 2 i = 10~ n s, r 32 = lCT 13 s and r 30 = 10 _12 s. 
The coupling constants in Eq. flUj) are a a = 10~ 4 C 2 /kg and = 5 • 10~ 6 C 2 /kg, and the 
linewidths of their corresponding transitions are T a = 2 ■ 10 13 s _1 and Ti, = I/V30 = 10 12 s _1 . 
Finally, the initial population density parameter is Nj nt = 5 ■ 10 23 m -3 . The population 
densities of the four level gain medium obeys the rate equations given in Eq. (jTTI) . and the 
gain medium interacts with the electromagnetic field through the gain polarization which 
obeys Eq. ( ITUl) . The Ag layer interacts with the electromagnetic field through the Drude 
polarization which evolves according to Eq. (J9j). 

Since the cylinder array is a single layer, it can be characterized as a metasurface- 18 . As 
a metasurface, the electromagnetic response is given by the surface polarizability 



Eq. ( l2Tj) is adapted from re£^, modified to be consistent with SI units and taking for granted 
that the metasurface is embedded in vacuum. The surface polarizability a is defined from 
the scattering matrix S. The S matrix is defined from the amplitude of the electric field of 
the scattered waves and is adjusted so that the effective thickness of the characterized array 
is zeroi 9 . For a symmetric and reciprocal array, such as the cylindrical spaser array, the S 
matrix components Sn = S22 are the reflection amplitude of a scattered wave and S12 = S21 
are the transmitted amplitude of the scattered wave. 

The surface polarizability of the cylindrical array is plotted in Fig. [2j The reflection 




(21) 



[1 + det(S) - (S n + S22)} e [(S 12 - S 21 ) - (S n - S 22 )} /c? 



(S12 - S21) + (Six ~ S22)} /c [1 + det(S) + (S u + S22)} ^ 
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FIG. 2: (a) Diagram of the simulation domain for the one dimensional cylindrical spaser array with 
a core gain medium (blue) and outer Ag shell (gray). A periodic boundary condition is imposed on 
the top and bottom boundaries, and a matched boundary condition (Sec. Ill C|) is imposed on the 
left and right boundaries. Real and imaginary parts of the electric surface polarizability a^ y (b) 
and magnetic surface polarizability a™" 1 (c) are plotted, clearly indicating separate electric and 
magnetic resonnaces. Inset are the field profiles for the two resonances and their corresponding 
wavelengths and Q factors. Color indicates magnetic field H^, and arrows indicate the electric 
polarization P = D — E. 

and transmission amplitudes used to calculate the surface polarizability were calculated 
from a frequency domain FEM simulation (COMSOL Multiphysics) where the Ag had a 
relative permittivity of €a 9 = 1 — u)p/(u>(u) — vy)) and the gain medium is simply a dielectric 
with permittivity = 9. We see from the surface polarizabilities that there is an electric 
resonance near Ao = 1220nm and a magnetic resonance near Ao = 830nm. Fig. [2] also show 
fields profiles for each of these resonances calculated using a FEM eigenfrequency simulation. 
Also shown are the wavelengths of the corresponding resonances A r = 27rc/Re(o; r ), and 
a resonance quality factor Q = 27rRe(u; r ) /Im(av) , where u r is a complex eigenfrequency 
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FIG. 3: (a) Electric surface polarizability a e y y and (b) magnetic sufrace polarizability a™" 1 for the 
cylindrical array with gain medium relative permittivity of = 9 — aTSSi n t/(u> 2 — u 2 — ir&u;). (c) 
Total, absorption as well as absorption in Ag and absorption in the gain medium. It is clear that 
the presence of the electronic transition in the gain medium strongly modifies the spectrum of the 
cylindrical array. 

returned by the same FEM eigenfrequency simulation. 

We are interested in using both resonances to achieve lasing, one resonance for enhancing 
the pumping of the gain medium and the other resonance for enhancing the lasing tran- 
sition. Therefore we choose the energy levels of the artificial four level gain medium so 
that the — ?- 3 transition approximately matches the higher frequency magnetic resonance 
Ub = 27rc/830nm, and the 1 — > 2 transition approximately matches the lower frequency 
electric resonance u a = 27rc/1221nm. The presence of a electronic transition will modify 
the spectrum of the cylindrical array for frequencies near that transition. Fig. [3] plots the 
surface polarizability near the magnetic resonance at Aq = 831nm for the cylindrical array 
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where the gain medium now has the relative permittivity ec = 9 — o&Nmt/ {u 2 — u\ — iT&u;). 
Fig. [3] also plots the total absorption of the cylindrical array, as well as separately plotting 
the absorption in the Ag and in the gain medium. Like Fig. [2J the data for these plots were 
calculated from a frequency domain FEM simulation. 

We can see from Fig. [3] that the interaction of the electronic transition with the magnetic 
shape resonance shown in Fig. [2j(c) causes these resonances to hybridize. As a result the 
response of the cylindrical array for frequencies near that transition is strongly modified. 
Instead of a single magnetic resonance we see now see multiple resonances, both electric 
and magnetic. Examining the absorption plotted in Fig. El^c) we see that the gain medium 
strongly absorbs at the magnetic resonance near Ao = 826nm. For our lasing simulations 
this will be the pump frequency. There is no way to know exactly what the lasing frequency 
will be without first running the time domain lasing simulation, except to say that it will 
be approximately equal to the frequency of the 1 — > 2 transition u a . A good initial guess is 
to set oji = u a , but after running the lasing simulation this can be adjusted to better match 
true lasing frequency. In what follows, we have used u± = 27rc/1219.3nm. 

Fig. H] shows data from a time domain simulation of the cylindrical spaser array using the 
parameters defined above. The initial state of the simulation is prepared with a previous 
simulation where the system is pumped with the field A2, with an intensity of 8W/mm 2 , 
while the incident probe field is set to A x = 0. The pump beam is turned on slowly with 
A2 having the profile 



where A pu is the amplitude of the pump beam, erf(x) = (2/a/F) J q x e~* dt is the error 
function, and r pu = 1.0 ■ 10~ 12 s is the pump rise time. The pump beam excites oscillators in 
the gain medium to the third energy level, which decays to the second energy level at the 
rate of l/r 32 . After t ^> r 2 i, the system is in steady state population inversion, but cannot 
lase since our model does not allow for generation of light due to spontaneous emission. The 
time is then reset, and the simulation shown in Fig. H] begins in this steady state population 
inversion. Shortly after t=0, a short probe pulse is emitted into the simulation domain. 
This excites the polarization field P^, which in turn begins the lasing process. The intensity 
of the resulting lasing field plotted in Fig. @Ja) spikes initially, but after about 30000 lasing 




(22) 
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t (2tt/ w i) x10 4 Pump Intensity (W/mm 2 ) 



FIG. 4: (a) Lasing intensity denned as the power emitted outward from the array in either direc- 
tion and (b) population inversion measured as J Q d 2 x(N y 2 — Nyi), where O is the domain of the 
simulation, and as a function of time normalized to lasing periods. The time domain simulation 
begins at t = with a steady state solution where the pump has been on for a very long time 
(t 3> T21) and the system has population inversion without lasing due to the lack of spontaneous 
emmision. (c) Plot of steady state lasing intensity vs. pump intensity. A linear fit indicates a 
pump threshold intensity of 7.15W/mm 2 and a slope of 0.145. 

periods it settles into steady state lasing. Fig. H] also plots the difference between the integral 
of the population densities N 2 and N x for the system beginning in population inversion. 

The time step used for the simulations in Fig. 4 varies throughout the simulation. When 
the pump is initially turned on the time step must be less then the pump rise time r pu . Once 
the pump is at a maximum the time step can be increased while the gain system approaches 
steady state. When the time is reset and a probe pulse is introduced the time step must 
be made smaller than the width of the probe pulse, and must remain small to resolve the 
resulting oscillations of the interaction between the probe pulse and the resonators as well 
as the initial exponential growth of the lasing beam. As the laser approaches steady state 
the time step can again be increased. At all times the time step must be smaller than the 
inverse rate of change of any transient beams (pump, probe or lasing beams). If oj\ is not 
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FIG. 5: Real (solid lines) and imaginary (dashed lines) parts of the electric surface polarizability af^ 
for the lasing electric resonance shown in Fig. [2] for various pump intensities (shown in legend). In 
Fig. 0[a) we see that at higher pump intensities the linewidth of the resonance narrows. In Fig. E^b) 
we see that at even higher pump intensities the imaginary part of the surface polarizability flips 
(indicating gain) and the linewidth of the resonance begins to broaden. 

close to the resulting lasing frequency the phase of Ai will rapidly change and will require a 
correspondingly small time step. Once the system begins lasing, the actual lasing frequency 
can be inferred from this oscillation in the phase of A 1; and U\ can be changed in the middle 
of the simulation. This causes the phase of Ai to slowly change allowing for a larger time 
step. 

There is a minimum pump intensity required for the light generated due to stimulated 
emission to overcome the internal losses in the cylindrical array. Fig. @Jc) plots the steady 
state lasing intensity vs. the pump intensity. A linear fit to the lasing data points indicates 
a threshold pump intensity of 7.15W/mm 2 . This threshold intensity depends on a number 
of variables, including all of the parameters of the gain medium, as well as the cylinder 
plasmon resonances used to enhance both the pump and lasing transition (Figs. [2] and [3]). 
These resonances in turn depend on the geometry and material parameters of the cylindrical 
array. 



17 



While there is a minimum threshold intensity for the array to exhibit lasing, we can 
observe interesting changes in the spectrum of the array at lower pump intensities. We saw 
by comparing Figs. [2] and E] that the spectrum of the cylindrical array was changed by the 
presence of the — > 3 transition. As we pump the array at increasing intensities we observe 
a similar change in the spectrum due to the 1 — > 2 transition. Fig. plots the surface 
polarizability (Eq. (|2T|) ) of the electric resonance for different pump intensities. These plots 
were created by pumping the cylindrical array with the field A2 for a long period of time 
(t ^> T21), and then injecting a Gaussian probe field Ai with a much weaker intensity. 
Applying a Fourier transform to the resulting time domain reflected and transmitted probe 
fields gives us the reflection and transmission amplitudes in the frequency domain, allowing 
us to calculate the surface polarizability according to Eq. ( l2Tj) . 

From Fig. [51 we see that for increasing values of the pump intensity, the lineshape of a™ 
resembles a Lorentzian 

ee ^0 fnn \ 

a W = °W - . .2 . ,2 ( 23 ) 

LO — 0J a — lJ a U) 

We see in Fig. [2(a) that as we increase the pump intensity, it is as if the positive valued 
scattering frequency 7 a is made smaller, narrowing the lineshape. We see in Fig. [2(b) 
that at even higher pump intensities, j a continues to shrink, passing through zero, and the 
imaginary part of changes sign, indicating gain. As the pump intensity continues to 
increase, 7 Q continues to grow more negative and the lineshape begins to broaden. 

Even though we have gain at the pump intensities in Fig. [3(b), we still do not have lasing 
because the gain is not large enough to overcome radiative losses. 



IV. CONCLUSION 



We have presented a finite element method simulation for a microphotonic lasing sys- 
tem. We have shown how to achieve a massive speedup in the simulation by separating 
the various fields into fields that oscillate at the carrier frequencies U\ or u; 2 , with slowly 
changing complex valued amplitudes. A demonstration of this simulation was provided with 
a two dimensional model of a one dimensional cylindrical spaser array as an example. The 
threshold pump intensity for this array was determined. Finally, we have shown how the 
linewidth of the lasing transition changes for various pump intensities. 
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